# download data https://africapolis.org
afripolis_shp = file.path(afripolis_dir,'africapolis')

# download data from ACLED
acled_fn <- paste0(wkdir,'/local/tab_data/Africa_1997-2020_Mar28-1.xlsx')

# download PETRODATA v11.zip data from PRIO.org https://www.prio.org/data/11  
# Lujala, P�ivi; Jan Ketil R�d & Nadia Thieme, 2007. 
# 'Fighting over Oil: Introducing A New Dataset', 
# Conflict Management and Peace Science 24(3), 239-256.
petro_shp = file.path(wkdir,'local/gis_data/016_society/PETRO_dataset_080907/PETRO_Onshore_080907')


require(Rcpp)
setwd(wkdir)
ntl_output_version <- '2020-08-04'
dist_output_version <- '2020-08-04'

calcPoints <- function(acled_loc){
  # add small value when no ntl #
  acled_loc[acled_loc==0] = 1
  rast.df <- as.data.frame(rasterToPoints(acled_loc))  
}


adm0_lake_chad_shp = paste0(wkdir,'/local/gis_data/003_boundaries/lake_chad/lake_chad_adm0')
adm0_reg_shp <- st_read(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
adm0_lake_chad_sf <- read_sf(dsn=dirname(adm0_lake_chad_shp), layer=basename(adm0_lake_chad_shp))
fishnet_ext <- extent(adm0_reg_shp)

# Maiduguri coordinates
maiduguri <- st_read(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/boko_haram_capital'),layer='Maiduguri_africapolis_pt')
bhstate <- st_read(dsn=paste0(wkdir,'/local/gis_data/003_boundaries'),layer='bhstate_ln')

## africapolis, poly
afripolis <- st_read(dsn=dirname(afripolis_dir),layer=basename(afripolis_dir))
afripolis <- subset(afripolis,ISO %in% c('CMR','NER','NGA','TCD'))
afripolis_pt = st_centroid(subset(afripolis,select=c("pop2015","pop1950","pop1960","pop1970","pop1980","pop1990","pop2000","pop2010")))

# select capitals from afripolis data #
capitals = c('Abuja','Niamey','Ndjamena/Kouss?ri [TCD]','Yaound?')
cap = subset(afripolis,Name %in% capitals)
cap_pt = st_centroid(cap)

# select primary cities from afripolis data #
pri = c('Lagos','Niamey','Ndjamena/Kouss?ri [TCD]','Yaound?')
pri = subset(afripolis,Name %in% pri)
pri_pt = st_centroid(pri)

# Petro
petro_ply = st_read(dsn=dirname(petro_shp), layer=basename(petro_shp),stringsAsFactors = FALSE)
oil_ply = subset(petro_ply,COUNTRY %in% c('Cameroon','Chad','Niger','Nigeria'))

##
## ACLED
##
acled.df <- read_excel(acled_fn)
print('subset for 4 countries')
conf.df <- subset(acled.df,COUNTRY %in% c('Nigeria','Cameroon','Chad','Niger'))

boko <- c('Islamic State (West Africa) and/or Boko Haram - Jamatu Ahli is-Sunnah lid-Dawatai wal-Jihad','Boko Haram - Jamatu Ahli is-Sunnah lid-Dawatai wal-Jihad')

conf.df <- as.data.frame(subset(conf.df, ACTOR1 %in% boko | ACTOR2 %in% boko))
conf.df$EVENTS <- 1

# gen combination indicator
conf.df$combconf <- conf.df$FATALITIES
conf.df$combconf[conf.df$combconf==0] <- 0.5

geoCRS <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
conf.spdf <- st_as_sf(conf.df, coords = c("LONGITUDE", "LATITUDE"), crs = geoCRS) 
conf.spdf$event <- 1
conf.spdf$EVENT_DATE = as.Date(conf.spdf$EVENT_DATE)
conf.spdf$month = as.integer(format(conf.spdf$EVENT_DATE,"%m"))
conf.spdf$year = as.integer(format(conf.spdf$EVENT_DATE,"%Y"))
conf.spdf <- subset(conf.spdf,year %in% 1999:2015)

xmin=0.1
xmax=24.1
ymin=1.6
ymax=23.6

# make an object the size and shape of the output you want
fishnet_bb <- matrix(c(xmin,  ymin,
                     xmax,  ymin,
                     xmax, ymax,
                     xmin, ymax,
                     xmin,  ymin), byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc(., crs = 4326)

# generate grid
grid01 <- st_make_grid(fishnet_bb, cellsize=0.1, crs = 4326, what = 'polygons', square = TRUE) %>%
            st_sf('geometry' = ., data.frame('ID' = 1:length(.)))
extent(grid01)


ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
if ((length(unique(ntl.zones$OBJECTID))==dim(ntl.zones)[1])==FALSE) { break }

##
## join pop data
##
require(dplyr)
pop.zones <- st_join(ntl.zones, afripolis_pt, join = st_contains)  %>% 
  group_by(OBJECTID) %>%
  summarize(pop1950 = sum(pop1950),pop1960=sum(pop1960),pop1970=sum(pop1970),pop1980=sum(pop1980),pop1990=sum(pop1990),pop2000=sum(pop2000),pop2010=sum(pop2010),pop2015=sum(pop2015))

# add area km2 #
ntl.zones$sf_area_km2 <- st_area(ntl.zones) / 1000000

##
## NTL DMSP
##
rprjlist=list.files(ntldmsp_dir, pattern=".tif$", recursive = TRUE, full.names=TRUE)

# load raster stack
ntl <- stack(rprjlist)
# set ntl na value 255 from .tif
NAvalue(ntl) <- 255
names(ntl) <- paste0("ntl_",substr(names(ntl),0,7))
geoCRS <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

# crop global to local region
ntl.crop <- crop(ntl,extent(grid01),snap='near')
extent(ntl.crop)
names(ntl.crop)
rm(ntl)

topcode=TRUE
if (topcode==TRUE) {
  ## add sum of top coded pixels ##
  print('top coded pixels')
  # reclassify the values into three groups 
  # all values >= 0 and <= 62 become 0, etc.
  m <- c(0, 62, 0,  62, 63, 1)
  rclmat <- matrix(m, ncol=3, byrow=TRUE)
  ntl.crop.topcode <- reclassify(ntl.crop, rclmat, include.lowest=FALSE)
  names(ntl.crop.topcode) <- paste0('ntl63sh',names(ntl.crop.topcode))
  ntl.topcode.extr <- exact_extract(ntl.crop.topcode,ntl.zones,fun=c('sum','count'),
                                    include_xy=FALSE, weights = NULL, progress=TRUE,max_cells_in_memory = 3e+7)
  ntl.topcode.dt = dplyr::bind_cols(ntl.zones,ntl.topcode.extr)
  ntl.topcode.dt$sf_area_km2 = st_area(ntl.zones) / 1000000
  ntl.topcode.df <- as.data.frame(ntl.topcode.dt)
  ntl.topcode.df <- subset(ntl.topcode.df,select=-geometry)
  names(ntl.topcode.df) <- gsub("\\.", "_", names(ntl.topcode.df))
  save.dta13(ntl.topcode.df, paste0(wkdir,"/proc_data/grid_ntltopcode",ntl_output_version,".dta"))
}


##
## continue to produce DIST.dta
##

##
## add area grid from raster
## 
ntl_area <- raster::area(ntl.crop[[1]])
names(ntl_area) <- 'r_area_km2'

ntl_stack <- stack(ntl.crop,ntl_area)

## stack
ntl.extr <- exact_extract(ntl_stack,ntl.zones,fun=c('mean','sum'),include_xy=FALSE, weights = NULL, progress=TRUE,max_cells_in_memory = 3e+7)
if (require(dplyr))
  ntl.dt = dplyr::bind_cols(ntl.zones,ntl.extr)

ntl.dt$sf_area_km2 = st_area(ntl.dt) / 1000000
ntl3.dt = ntl.dt

##
## add distance to poly to center point
##

albers_africa = "+proj=aea +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"
ntl.zones.prj <- ntl.zones %>%
  # 102022 = WGS 1984 Albers for Africa
  st_transform(albers_africa) %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)

## boko haram area/polygon
bhbr = read_sf(dsn=paste0(wkdir,"/local/gis_data/003_boundaries"),layer="bhbr_dst")

## dist from nngeo ##
(require(nngeo))
ntl3.dt$dist2bordbhbr <- unlist(nngeo::st_nn(ntl.zones.pts, bhbr, k = 1, returnDist = T)$dist) / 1000
ntl3.dt$dist2cap_nga <- unlist(nngeo::st_nn(ntl.zones.pts, cap_pt[cap_pt$ISO=='NGA',], k = 1, returnDist = T)$dist) / 1000
ntl3.dt$dist2pri_nga <- unlist(nngeo::st_nn(ntl.zones.pts, pri_pt[pri_pt$ISO=='NGA',], k = 1, returnDist = T)$dist) / 1000

dist.dt <- as.data.frame(ntl3.dt)
dist.dt2 <- dist.dt %>%
  dplyr::select("OBJECTID",
                "dist2cap_nga",
                "dist2pri_nga",
                "dist2bordbhbr")

save.dta13(dist.dt2, paste0(wkdir,"/proc_data/dist",dist_output_version,".dta"))